Interplay of local hydrogen-bonding and long-ranged dipolar forces in simulations of confined water 
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Spherical truncations of Coulomb interactions in standard models for water permit efficient molecular sim- 
ulations and can give remarkably accurate results for the structure of the uniform liquid. However truncations 
are known to produce significant errors in nonuniform systems, particularly for electrostatic properties. Local 
molecular field (LMF) theory corrects such truncations by use of an effective or restructured electrostatic poten- 
tial that accounts for effects of the remaining long-ranged interactions through a density-weighted mean field 
average and satisfies a modified Poisson's equation defined with a Gaussian-smoothed charge density. We apply 
LMF theory to three simple molecular systems that exhibit different aspects of the failure of a naive applica- 
tion of spherical truncations - water confined between hydrophobic walls, water confined between atomically- 
corrugated hydrophilic walls, and water confined between hydrophobic walls with an applied electric field. 
Spherical truncations of 1/r fail spectacularly for the final system in particular, and LMF theory corrects the 
failings for all three. Further, LMF theory provides a more intuitive way to understand the balance between 
local hydrogen bonding and longer-ranged electrostatics in molecular simulations involving water. 



I. INTRODUCTION 

An accurate and efficient treatment of Coulomb interactions 
in molecular simulations represents an important and ongoing 
challenge. Standard biomolecular simulation packages like 
CHARMM jUl and AMBER |2l quite generally assign effective 
point charges to interaction sites even in neutral molecules 
in order to approximate the charge separation of the elec- 
tron cloud along polar bonds. Effective point charges are also 
found in most standard water models such as the extended 
simple point charge (SPC/E) model 13] shown in Fig. la, and 
water molecules are increasingly being included explicitly in 
biomolecular simulations. 
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FIG. 1: Bulk SPC/E water at 300 K and 0.998 g/cm^. (a) Diagram 
of SPC/E model of water. Point charges are assigned to each atomic 
site and a Lennard- Jones potential is centered about the oxygen with 
core diameter au indicated by dashed lines, (b) Various site-site pair 
correlation functions in bulk SPC/E water are each displaced verti- 
cally by 0.5 for clarity. Gaussian-truncated (GT) water using vo{r) 
with (T = 0.40 nm agrees well with a full treatment of electrostatics, 
(c) Deviations (5gQQ, defined as the difference between g(r) for GT 
water and g(r) using full electrostatics. Simulation error bars are also 
shown. 



In order to minimize edge effects in necessarily finite sim- 
ulation cells, these simulations use periodic boundary condi- 
tions 101. Researchers have long sought to devise spherical 
truncations of the Coulomb 1/r potential so that the truncated 
potential accurately describes the strong Coulomb forces be- 
tween closely-spaced charges that lead to hydrogen bonding 
in water but then vanishes quickly beyond some properly cho- 
sen cutoff radius. Then, only the minimum (closest) image of 
an interacting charge needs to be accounted for and fast and 
efficient simulations that scale linearly with system size are 
possible. In bulk-like systems these spherical truncation ap- 
proaches can be surprisingly accurate, but standard estimates 
of thermodynamic properties in the bulk liquid are less satis- 
factory illiaili[ii[ni[il. 

However, it has long been established ifisll that when these 
spherical truncations are employed in nonuniform geometries, 
such as near a Upid bilayer, there can be pronounced er- 
rors in structure, thermodynamics, and particularly electro- 
static properties. In those cases, the neglected long-ranged 
forces tl4i] compete with the local ordering due to cores and 
hydrogen bonding captured well by spherical truncations. 

These failures have led most researchers to conclude that 
the success of spherical truncations in bulk was seductive but 
misleading. Instead, lattice summation techniques - the elec- 
trostatic method of choice - give an accurate account of all pe- 
riodic images. Ewald [15] and Lekner lfl6ll summations have 
long been used in simulations of charged and dipolar fluids 
and biomolecules. Much current research is focused on devel- 
oping more efficient implementations such as particle-mesh 
Ewald |[Il[ll and in removing some of the artifacts that may 
arise from the use of periodic lattice sums in fluid environ- 
ments lfl9ll . 

But is it possible to build on the success of spherical trun- 
cations of 1/r interactions in uniform environments lIH [tI 
H H [H [nl [H and to understand and correct their fail- 
ures in nonuniform environments |[l3l We show that 
local molecular field (LMF) theoryloll IH [H offers 
exactly that opportunity. LMF theory corrects naive trun- 
cation schemes by introducing an effective electrostatic po- 
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tential Vn{r) that accounts for the remaining effects of the 
long -ranged interactions. LMF theory also suggests a partic- 
ularly effective smoothed Gaussian-truncation scheme for the 
charges in SPC/E water and other interaction-site molecular 
models. Furthermore, the theory provides a general concep- 
tual framework that permits a qualitative as well as quanti- 
tative understanding of the effects of long-ranged Coulomb 
interactions on local water structure - something that is often 
obscured in black box application of lattice summation tech- 
niques. 

We will briefly discuss the application of the Gaussian- 
truncation scheme to bulk SPC/E water. Then we will demon- 
strate that the inclusion of the effective potential V/j(r) given 
by LMF theory is the sole modification of molecular dynam- 
ics simulations of Gaussian-truncated (GT) water required 
to yield accurate minimum image simulations of water in 
nonuniform slab environments. While we focus on systems 
composed of solely water molecules and simple confinement 
potentials, we argue below that the LMF approach is relevant 
for the treatment of electrostatics in general biomolecular sim- 
ulations. 



II. BULK GAUSSIAN-TRUNCATED (GT) WATER 

In LMF theory, we split the 1/r potential associated with 
each charge into rapidly-varying (vq) and slowly-varying (vi) 
components as 



1 crfc(r/(T) crffr/a) , , , , 

- = '-^ + =Vo{r)+vi{r). (1) 



This potential separation isolates strong short-ranged and 
rapidly-varying interactions in vo{r), and the remaining 
slowly-varying long -ranged forces are encompassed by vi (r). 
This interaction vi (r) is proportional to the electrostatic po- 
tential arising from a smooth normalized Gaussian charge dis- 
tribution with width <t, and is defined via the convolution 
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Thus vi (r) is slowly-varying in r-space over the smoothing 
length a and is simultaneously localized to small wave vectors 
in reciprocal space, as can be seen from its Fourier transform 
vi{k) = 47rfc-2cxp[-(fccr)V4]. 

Hence vo{r) = 1/r — vi (r) is the screened potential result- 
ing from a point charge surrounded by a neutralizing Gaussian 
charge distribution whose width a also sets the scale for the 
smooth truncation of vq. For small r < a, the force due to 
vo{r) approaches that of the bare point charge. By increasing 
a we increase the effective range of essentially unscreened 
Coulomb interactions included in vo{r). 

LMF theory introduces a general mapping from the original 
system with full Coulomb interactions and an applied electro- 
static potential V(r) arising from fixed charges or an applied 
electric field to a well-chosen "mimic system" with Coulomb 



interactions replaced by the short-ranged vo{r): 
1/r 1 r Voir) 
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Here, Vii{r) is a restructured electrostatic potential in the 
mimic system that accounts for the averaged effects of the ne- 
glected long-ranged forces represented by in essence by 
solving a modified Poisson's equation as explained in a later 
section. 

In this paper, we simulate the SPC/E model of water as 
shown in Fig. la. All data are results from molecular dy- 
namics simulations as described in the Methods section. The 
point charges in SPC/E water play a dual role, generating both 
strong short-ranged forces that lead to the local hydrogen- 
bond network as well as the longer-ranged dipole-dipole in- 
teractions. Neither the point charges nor the Lennard-Jones 
(LJ) potential describing the molecular cores are perfect rep- 
resentations of actual water interactions but adjusting the con- 
stants makes the total potential quite reasonable for tempera- 
tures and densities of interest. 

Since the charges generating the local hydrogen bonding 
network in water are very strongly coupled, a reasonable first 
approximation is simply to replace all point-charge interac- 
tions by the truncated vo{r) without employing a Vr. LMF 
theory suggests that Coulomb interactions are best truncated 
on a site-site basis, since it is the basic 1/r potential that we 
attenuate about each site. Such a site-site scheme has been 
implemented in other spherical truncations of 1/r interactions 
as well f lO']. When the wo(r) truncations are made, we have 
Gaussian-truncated (GT) water. 

Crucial to the success of the Gaussian-truncation scheme 
is the choice of a cr that will be large enough to accurately 
describe the nearest-neighbor hydrogen bonds. Using a cr of 
0.4 nm or larger, GT water can give a remarkably accurate 
description of all O and H pair correlation functions in bulk 
water, as shown in Fig. lb. The difference between gQQ from 
a full electrostatic treatment (indicated by the labels "Full" in 
the figures) and ggo from GT water is shown in Fig. Ic to 
emphasize the good agreement. We expect that certain earlier 
truncation schemes iflOll with very slowly-varying vi{r) could 
give comparably accurate results. 

The local hydrogen-bond network is often viewed as the 
most important qualitative feature of water ll24ll . and GT wa- 
ter captures that feature with a very good description of near- 
est neighbor hydrogen-bond energetics and local tetrahedral 
order, while ignoring the longer-ranged electrostatic effects of 
the bound charges. As argued elsewhere ||25||, the accuracy of 
GT water in bulk arises from a strong cancellation of the long- 
ranged electrostatic forces exerted by surrounding molecules 
in the uniform bulk on any given molecular site. These site- 
based spherical truncations can also yield highly accurate an- 
gular correlation functions ll25ll . a property that is missed with 
molecule-based truncations that do not neglect only slowly- 
varying forces i^. LMF theory also suggests simple analytic 
corrections to the energy and pressure such that the thermo- 
dynamics of bulk GT water are accurate as well ll26ll . 
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FIG. 2: Model hydrophobic wall confinement of SPC/E water with 
an applied field. GT water with cr = 0.6 nm fails while treatment 
with full LMF theory succeeds, (a) The simulation system and near- 
est images in the y-direction confined by the Leimard- Jones wall po- 
tential, (b) The oxygen density profile relative to the bulk density, 
(c) The charge density profile. Density profiles are only shown for 
z < —1.0 nm for greater clarity of the atomic level behavior, though 
these densities are strictly asymmetric about the origin. 



LMF-CORRECTED GT WATER IN AN APPLIED 
FIELD 



m. 



Following previous work 11271 12811 . we confine SPC/E wa- 
ter between smoothed hydrophobic LJ walls and apply a con- 
stant electric field Eq normal to the walls, as sketched in 
Fig 2a. The confined water should behave like a dielectric 
slab from introductory electrostatics, resulting in a polariza- 
tion field Epoi opposing the applied electric field 1221 ■ The 
oxygen density profile and charge density profile using only 
GT water are completely incorrect when we apply an electric 
field of 10.0 V/nm as shown in Fig. 2b and Fig. 2c. 

Using spherical truncations alone neglects the long-ranged 
electrostatic consequences of the local ordering of the wa- 
ter molecules, which act to attenuate and weaken the applied 
electric field in the central region as suggested by Le Chate- 
lier's principle. Thus the GT water molecules spectacularly 
over-order with the unattenuated electric field since the ex- 
ternally applied field induces an incorrect, strongly polarized 
state throughout the GT water slab that happens to be compat- 
ible with an energetically-favored arrangement of hydrogen 
bonding. Inclusion of a self-consistent Vu{z) to represent the 
long-ranged net effect of water molecule polarization makes 
the density profiles again agree very well with the results for 
the full system as shown in Fig. 2b and Fig. 2c. With the 
inclusion of net long-ranged effects, the highly-ordered state 
does not occur because it is entropically disfavored and thus 
is highly unlikely in the true attenuated field. 

The inclusion of a restructured Vij(z) also corrects many 



other properties of confined water when an electric field is ap- 
plied such as the dielectric constants and details of molecular 
orientation profiles ll26ll . However, the density profile is the 
simplest and most spectacular example of the failure of GT 
water in this system, and LMF theory handily corrects that 
failure. 



IV. LOCAL MOLECULAR FIELD THEORY FOR 
ELECTROSTATICS 



LMF theory has been applied to many different ionic sys- 
tems, including charge mixtures ll23l Isol IsTl [32ll . Gen- 
eral derivations of the approach are available in the litera- 
ture 1I20I I21I I22I1 . However, when the LMF truncation is 
applied to the general 1/r functional form and the same a 
is applied to every charge-charge interaction in the system 
regardless of identity, we may write a substantially simpler 
formulation of the LMF mixture equation than previously 
stated ifsoll . An explanation of these symmetries and the va- 
lidity of the LMF equation below for both charge mixtures 
and site-site molecular models will be detailed in forthcoming 
papers iHH. 

Provided that the cr chosen is at least large enough for 
vo{r) to incorporate the nearest-neighbor interactions such as 
hydrogen-bonding in water and for vi (r) to be slowly varying 
over those nearby correlations, LMF theory predicts that the 
restructured potential Vr is accurately given by 



Vnir) = V(r) 
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The notation p'j^{r; [Vb] ) simply indicates that we now use the 
charge density profile of the mimic system (indicated by R) 
in the presence of the restructured potential Vr. The charge 
density p'^ encompasses the effect of vo{r) and Vr as well 
as all other non-electrostatic interactions like the hydrophobic 
walls and the LJ cores in SPC/E water. We use the constant 
C to set the zero of our potential. In the systems examined in 
this paper, this formula may be rewritten as a one-dimensional 
equation as detailed in Ref. ll30ll . 

Vfl can be understood as a mean-field electrostatic poten- 
tial representing the long-ranged forces that may be reason- 
ably averaged over. Alternately, harnessing the convolution 
definition of vi{r) in Eq. |2l we may rewrite the LMF equa- 
tion as the convolution of a smoothed charge density with the 
full Coulomb potential. Thus, Vn exactly satisfies a modified 
Poisson's equation 

{Vfl(r) - V(r)} = --pl'ir; [V^]), (5) 
Co 

defined with a Gaussian-smoothed charge density 
1 
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This again emphasizes that all charge interactions in the sys- 
tem must be identically truncated for this modified electrostat- 
ics to be valid. 
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The form of the short-ranged potential vo{r) might suggest 
an analogy to Ewald summation. However, LMF theory fo- 
cuses on the short-ranged mimic system and how best to use 
and correct its properties based on results from equilibrium 
theory, while Ewald summation focuses on the full system 
with exact periodic boundary conditions. As a result, is 
a static external potential based on equilibrium properties that 
is determined self-consistently as described in the Methods 
below, whereas Ewald summation evaluates fluctuating forces 
arising from both the short- and long-ranged components of 
the Coulomb interactions at each time step. 



V. CLASSICAL FAILURE OF GT WATER EXPLAINED 

While we believe that the above failure of GT water alone in 
an electric field has not been noted before, problems with elec- 
trostatic properties of confined truncated water in zero field 
are well known lisll . Shown in Fig. 3 is the classic example 
of the failure of GT water in nonuniform situations - the in- 
correct induced potential profile of water molecules. SPC/E 
water is again confined between the hydrophobic walls, now 
with no applied electric field. As well-established by Rossky 
and coworkers ^T^, water molecules next to the wall tend to 
point one hydroxyl group into the wall so that one hydrogen 
bond per molecule on average is broken rather than two. This 
induces a dipole layer near each hydrophobic surface with 
the dipole pointing toward the walls. This atomically-detailed 
dipole layer leads to a semi-discrete jump at the surface in the 
electrostatic potential defined as 

$poi(^) = -- r r dz"p'^{z"). (7) 

£0 J-L/2 J-L/2 

This potential $poi should plateau in the central bulk region; 
however, the use of GT water fails dramatically in this re- 
spect. Due to the initial overorientation of water molecules at 
the surface, driven by the maintenance of hydrogen bonds in 
GT water, the polarization potential never plateaus. This fail- 
ure is well documented ifTsl [Till and is attributable to the ne- 
glect of net-additive long-ranged forces due to the formation 
of dipole layers at the surface. As also shown in Fig. 3, LMF 
theory correctly includes those neglected forces with tremen- 
dous success. 

The form of Vji{z) shown in Fig. 4 highlights its physi- 
cal effect on GT water molecules. Without Vr, the truncated 
molecules overorient near the wall. The restructured electro- 
static potential Vr from LMF theory reflects the net long- 
ranged forces due to the ordering of water near the surfaces 
by applying a reorienting torque on surface water molecules. 
Since Vr is defined via a self-consistent equation, it repre- 
sents a statistical mechanical balance between the favorabil- 
ity of maintaining local hydrogen bonds and the electrostatic 
penalty for creating an overly severe dipole layer The po- 
tential Vr varies smoothly with z because it is determined by 
averaging over the slowly-varying vi{r) in Eq. |4] This sug- 
gests that a simple theory yielding this general form may be 
possible. 




FIG. 3: Polarization potential profile,<I>poi(2;), between hydrophobic 
walls with no applied field. GT water does not display the expected 
plateau in the bulk region. Treatment of the system with LMF theory 
yields strong agreement with full treatment of the electrostatics. 
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FIG. 4: FoiTn of Vr{z) for water confined between hydrophobic 
walls with no applied electric field. The self-consistent potential ap- 
plies a reorienting torque on water molecules near the surface, mim- 
icking the effect of the long-ranged dipole layers. 



While the physics of this failure was convincingly estab- 
lished by Spohr its structural roots in the charge den- 
sity profile remained mysterious. Brooks and coworkers ob- 
served an overorientation of molecular dipoles near the walls 
in spherically truncated water lfl3ll . But as shown in Fig. 5a, 
the charge density profile for GT water tracks quite closely 
with the results from full Ewald summation techniques for 
SPC/E water A dipole layer near the wall seems quite clear, 
and there are only slight deviations in the peaks next to the 
walls. However, the polarization potential $poi belies the ap- 
parent accuracy of the atomic-level charge density profiles for 
GT water. The polarization potential may be fit accurately by 
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FIG. 5: Charge density profiles for water confined between hy- 
drophobic walls with no applied field, (a) The atomic level charge 
density for 2 < 1.0 nm shows only slight discrepancies between 
GT water and either full or LMF-treated water. This density profile is 
symmetric about the origin, (b) The gaussian-smoothed charge den- 
sity p'" shows the failure of GT water to capture the true electrostatic 
dipole layer much more clearly. In either instance, the LMF-treated 
water is in very good agreement with a full treatment of electrostat- 
ics. 



a quadratic function in the central region between the walls. 
This fit suggests the existence of a constant residual charge 
density of approximately —0.04 Co/nm"^ throughout the "bulk" 
region neutralizing the overaccumulation of positive charge 
next to each wall. 

LMF theory allows us to highlight exactly that relevant 
structural problem. The writing of the LMF equation as a 
modified Poisson's equation in Eq.|5] suggests that the charge 
density that best reflects the long-ranged electrostatics is not 
p'i, but rather the Gaussian-smoothed charge density p'^" . The 
smoothed charge density p'*" in Fig. 5h clearly shows the 
dipolar layer at each wall for the water confined with no 
electric field. Beyond that, GT water alone very obviously 
spreads the dipole layer into the bulk. The residual nega- 
tive charge density is approximately the value suggested by 
the form of $poi for GT water. Since the smoothed charge 



density is two orders of magnitude smaller than the atomic 
level charge density, this underlying electrostatics is buried 
alongside the slight simulation noise in p'^. On Gaussian- 
smoothing, the random simulation noise in p'^ cancels, un- 
covering the smaller but coherent electrostatic features. In 
a sense, the strong charge layering over 0.2 nm displayed in 
Fig. 5fl is not the long-ranged dipole near the wall. The true 
dipole layer is due to average behavior over a few molecular 
layers. Classical electrostatic textbooks often posit some sort 
of smoothing over such local order 1I29I1 . 

This smoothed charge density allows for a more direct anal- 
ogy between molecular simulations and classical electrostat- 
ics. The LMF approach naturally distinguishes the long- 
ranged charge-charge interactions representing overall dipolar 
ordering from the short-ranged charge-charge interactions in 
simulations that capture hydrogen-bonding and short-ranged 
polar attraction. Given that any g greater than some CTmin 
will describe relevant short-ranged interactions like hydrogen 
bonding and can be used in simulations, too much significance 
should not be assigned to the particular chosen <t value or the 
exact magnitudes of the resulting p^" . However we strengthen 
a previous suggestion 13411 that charge density profiles reveal 
electrostatic effects in a less-biased fashion than molecule- 
based profiles: perhaps a smoothed profile is even less biased. 



VI. LMF THEORY IS MORE BROADLY APPLICABLE 

We now illustrate the application of LMF theory to mod- 
erately more realistic molecular surfaces by replacing the hy- 
drophobic walls with empirical Pt(l 1 1) surfaces that order the 
water molecules at the surface, attracting the oxygen atoms to 
localized binding sites as shown by Berkowitz and cowork- 
ers llssll . This surface is meant to represent specific interac- 
tions and detailed ordering of water molecules at an atomi- 
cally corrugated surface, thus mimicking a hydrophilic inter- 
action. The particular wall potential we use does not have ex- 
plicitly fixed wall charges, but LMF theory can readily handle 
these as well lf33ll . 

Since this surface induces a charge density that is explicitly 
a function of x, y, and z, one might reasonably conclude that 
we now must solve the full three-dimensional LMF equation 
given in Eq. |4] to obtain a Vr that is a function of r. How- 
ever, the charge smoothing inherent in the LMF equation jus- 
tifies a much simpler treatment of Vr. The lateral variations in 
charge density are washed out when smoothed with the Gaus- 
sian of width cr, leaving only a slowly-varying p"^" (z) in the 
Poisson-Iike formulation of the LMF equation given in Eq. |5] 
This p"^" (2) again reflects the formation of dipole layers near 
each wall. Thus, solving a one-dimensional LMF equation 
yields remarkably accurate results for the polarization poten- 
tial profile <l>poi(2), as shown in Fig. 6. Such findings give 
us confidence that the LMF treatment of interesting and more 
complex simulations such as water-membrane interfaces and 
liquid-liquid interfaces will prove to be just as simple. 
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FIG. 6: Polarization potential profile, $poi(z), between corrugated 
walls with no applied field. GT water again does not display the 
expected plateau in the bulk region. Treatment of the system with 
LMF theory again yields strong agreement with full treatment of the 
electrostatics. 



VII. CONCLUSIONS AND OUTLOOK 

Here we have demonstrated the utility of the LMF approach 
both in simulating a molecular system in a slab geometry and 
also in analyzing the structural properties of this system. The 
inclusion of a restructured electrostatic potential Vr allows for 
remarkably accurate results in this slab geometry. The spec- 
tacular failure of the spherical truncation of 1 /r interactions 
to describe even the density profile of water in an applied field 
is corrected through the use of LMF theory. Earlier known 
failures of truncated water models to describe the electrostatic 
potential for confined water in zero field are also corrected 
by LMF theory regardless of the specific details of the con- 
fining surfaces. Furthermore, the splitting of 1/r into short- 
ranged and long-ranged components provides a natural way 
to disentangle the impact of local hydrogen-bonding and of 
long-ranged dipolar interactions in various aqueous systems. 
In light of LMF theory's success in this slab geometry, we 
certainly believe that it could be of use in other physically 
relevant slab-geometry simulations such as biological mem- 
branes and liquid-liquid interfaces. Perhaps it also could help 
resolve a recent dispute over whether or not depletion of wa- 
ter between capacitor plates occurs in the simulation of an 
open system llSq, [37i [ssll . Furthermore, given previous work 
treating primitive models of ionic solutions with LMF the- 
ory 1231 [30l [32I1 . the study of aqueous ionic solutions is an- 
other natural step forward. Truncating all 1/r interactions via 
vo{r) with the same a would lead to a natural mimic system 
for such ionic solutions. 

Examination of either Vr or the smoothed charge density 
yo'" prescribed by LMF theory also clarifies the underlying 
electrostatics and the failure of spherical truncations. The 
Gaussian-smoothed perspective for long-ranged electrostatics 



greatly elucidates the perceived contradiction between seem- 
ingly good charge density profiles but very poor polarization 
potential profiles for water confined between walls in zero 
field. Analysis of long-ranged electrostatic properties using 
the smoothed charge density seems more plausible because 
the point charges in classical site-site models approximate the 
real quantum-mechanical electron cloud, and longer-ranged 
electrostatic effects should be relatively insensitive to the ex- 
act placement and magnitude of point charges in models with 
the same dipole moment. Thus in GT water, the short-ranged 
cores are altered to reflect both hydrogen-bonding and polar 
attractions represented by vo{r) as well as the LJ cores. Sim- 
ulation using GT water without a self-consistent will only 
be reasonable if the resulting p''" reflects the underlying elec- 
trostatics. Since Vr is the solution to a modified Poisson's 
equation, we should not expect this to be generally true. 

Overall the use of LMF theory for electrostatics in sim- 
ulations seems very promising, but certainly this approach 
needs more investigation. Upcoming articles will detail the 
derivation of the simple LMF electrostatic potential given in 
Eq.|4|for charge mixtures Issll and for site-site molecular mod- 
els I25I as well as simple anaWtical thermodynamic coiTec- 
tions for these mimic systems ll26ll . Further investigation into 
optimizing the solution of the LMF equation is called for, 
and the slowly-varying nature of the smoothed charge den- 
sity suggests that simple approximations could prove quite 
useful I32il . At present, the LMF mapping is for equilibrium 
properties only. However current work has shown that LMF 
theory used in tandem with with a Born-Oppenheimer-like ap- 
proximation can describe the slow dynamics and the associ- 
ated free energy landscape of charged polymers folding in the 
presence of mobile salt ions lf39ll . In general, however, a full 
extension of these ideas to dynamics remains a challenging 
open question. Its resolution may clarify the relationship be- 
tween the fluctuations of these short-ranged systems and those 
in experimental systems. 

While all nonuniform systems examined herein have a slab 
geometry, LMF theory is valid for general nonuniformity. As 
such, we view the work presented in this paper as a demon- 
stration of the utility of LMF theory to treat and to understand 
electrostatic interactions in molecular systems far more com- 
plex than those studied here. 



VIII. METHODS 

For this research, we have used the dlpoly2. 16 molecular 
dynamics package ll40ll modified to include the Lennard- Jones 
wall potential llj, the Pt(lll) wall potential js?], the LMF 
theory vo{r) potential and the Vr{z) mean field potential, as 
well as the three-dimensional corrected Ewald sum for slab 
systems (EW3DC) iH. SPC/E water il) was simulated at 
constant volume and temperature using the Berendsen ther- 
mostat |42| with a timestep of 1 fs. For the bulk water pair 
correlation functions, we simulated 1728 SPC/E molecules at 
bulk density for 1.5 ns after 500 ps of equilibration at a tem- 
perature of 300 K. For each nonuniform system and electro- 
static technique we also have carried out at least 500 ps of 
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equilibration and subsequently 1 .5 ns of simulation at 298 K. 
The lateral simulation box size for the smooth Lennard-Jones 
walls was 2.772 nm x 2.772 nm with 1024 SPC/E molecules, 
and the walls were placed at ±2.25 nm in order to obtain the 
bulk density of water in the central region. When the elec- 
tric field of 10.0 V/nm was applied, the walls were moved 
to 2.153 nm to maintain the bulk density of water. (There- 
fore effects of electrostriction were not accounted for) The 
Pt(lll) wall area was 3.05 nm x 2.88 nm with 1054 SPC/E 
molecules, and the for the walls (as defined in Ref. Issll ) 
were placed at ±1.56 nm. In both cases the simulation cell 
spanned 14.0 nm in the z-direction, though this was only cru- 
cial for the Ewald summations to maintain sufficient vacuum 
space between periodically replicated slabs. 

The three electrostatic techniques used for the slab geome- 
try were Ewald summation ifisll . minimum image simulations 
using a spherically-truncated vq (r), and minimum image sim- 
ulations using a spherically-truncated vo{r) and Vfi.{z)- The 
Ewald technique used was the corrected three-dimensional 
Ewald scheme of Yeh and Berkowitz ll4lll . We used an a of 
3.40 nm^^ and fc- vectors spanning (12, 12,60) for LJ walls and 
( 1 3, 12,60) for Pt( 1 1 1) walls. For vq (r) we used a <t of 0.60 nm 
for slab simulations and 0.40 nm for the bulk simulations. For 
bulk, using a larger <t of 0.60 nm produced equivalent results, 
though the 5g will be slightly different. 

The LMF equation must be solved self-consistently so that 
the determined via simulation using Vfl yields the same 
LMF equation prediction of Vn- Self-consist ently solving the 
LMF equation in previous simulation work pill proved rel- 
atively straightforward with a linear mixing parameter A as 
Vi+^ = A ■ V^^^ + (1 - A) • V), where V)f ^ is the restruc- 
tured potential determined by inserting the p'^ due to V]^"^ into 
Eq. |4] However, for these systems, the solution proved more 
challenging due to the need to distinguish between inherent 
equilibrium fluctuations in molecular orientations in the short- 
ranged system and the induced response of the short-ranged 
system to a variation in the small and slowly-varying Vr. 

We may still solve the LMF equation fairly quickly, given a 



reasonable initial Vj^., by carrying out A'^ simultaneous short- 
time simulations with distinct starting points during the evo- 
lution to a self-consistent solution. The averaged p"^ of the N 
systems simulated in V|j at the i*'* solution step better revealed 
the response of the system to the change in Vr from V]^^ to 
V]j. These short-time simulations each ran for 25 ps of equi- 
libration and 50 ps of accumulation with a larger At of 2.5 fs. 
Then the p'^{z) resulting from each of the N short-time sim- 
ulations are averaged together and inserted in the LMF equa- 
tion to yield V]|^^. Just as for the simpler solution schemes, 
we then used linear mixing with A ranging from 0.05 to 0.10 
to predict the subsequent V]^^. For a ~ 0.60 nm, Vr con- 
verged within 3 steps of 10 parallel simulations with no ap- 
plied Eq and within 10 steps of 10 parallel simulations with 
an applied field of 10.0 V/nm. However, this process was em- 
ployed simply to demonstrate the utility of a self-consistent 
LMF treatment of electrostatics in molecular simulations. An 
optimized approach will likely appear quite different. Cur- 
rently, using the interpretation of the LMF equation given in 
Eq. |5] we are exploring options such as fast Poisson solvers 
and Car-Parinello techniques in solving the equation adap- 
tively during simulation. Such solutions for would still 
be static in nature once self-consistency was achieved. 
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